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Abstract. - We use simulations to investigate collision time distributions as one approaches 
the static limit of steady-state flow of dry granular matter. The collision times fall in a power- 
law distribution with an exponent dictated by whether the grains are ordered or disordered. 
Remarkably, the exponents have almost no dependence on dimension. We are also able to 
resolve a disagreement between simulation and experiments on the exponent of the collision 
time power-law distribution. 



Dense granular matter does not easily fit into our standard classification of matter. Interest 
in the physics community was stimulated by the observation that forces between particles in 
the system are exponentially distributed suggesting that any load on the system is carried by 
a small number of force chains [1,2]. More recent studies have suggested that spatial ordering 
is a key factor in the force response [3,4], something that may not be clearly distinguishable 
from the exponential tail of the force distribution. Further work, both experimental and 
theoretical, has suggested that perhaps it is actually the low end of the force distribution that 
should be examined when deciding whether a system is jammed [5-7]. 

As the distribution of forces in a static granular system is history dependent [8] , it makes 
sense to explore the static limit of dynamic models in order to shed light on the origin of 
jamming. In this letter, we perform simulations of gravity driven dense granular flow in 
two and three dimensions using an event-driven model involving only binary collisions. We 
reproduce the known power law relationship between the mean flow velocity and velocity 
fluctuations [9] and explain the discrepancy between the experimentally observed power law 
for the distribution of collision times [10] and that found in previous simulations [11]. Further, 
we find that the distribution of collision times is capable of differentiating between ordered 
and disordered glassy systems. Moreover, the exponents appear to be superuniversal in the 
sense that they are independent of dimension. In addition, we find that the force distribution 
turns up at low impulses in the jammed state. 

We simulate hard spheres as they fall down a rectangular chute under the influence of 
gravity as shown in Fig. ^a). At the bottom of the chute a sieve is modeled by having 
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Fig. 1 - (a) Section of a simulation involving 43 200 grains with 1 % polydispersity in a 32a x 32a x 400a 
system. Average (b) density for 1% (dashed line) and 15% polydispersity (solid line) along the height 
of a 3D chute. For 15% polydispersity, (c) the y— velocity and (d) average acceleration, in units of 
g and as measured by the material derivative dv y /dt — dtv y + v a d a v y . The inset in (d) shows the 
absolute value of the same data on a log scale, showing the lg acceleration in the free fall region, (e) 
shows the pressure tensor components P xx as the solid line and P yy as the dashed line. The inset in 
(d) shows that the pressure approaches Eq. (|3J as jamming is approached from the fluid phase [4]. 



the particles reflecting from the bottom with a probability p (typically p — 10%). Particles 
transmited through the bottom reappear at the top of the chute once again to fall down 
through it. Particles reflect off the walls of the chute with a partial loss, typically 10 %, 
in their vertical velocity. This is a simple effective way to model rough walls. The grain 
polydispersity is varied from to 15% in different simulations [12]. 

Velocities after collision ri and r 2 in terms of velocities before collision, ri and r 2 , are 



( ri \ _ / h V (1 + M) f -mi m 2 \ f f r q \ 
\ ^2 J V *a / (mi+m 2 ) V TO i _m i / V r 2-q / 



(1) 
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where q = (r 2 — ri)/|r 2 — ri| [13,14]. p, is a velocity-dependent restitution coefficient described 
by the phenomenological relation [15, 16], 

/ \ _ | 1 - (1 - Mo) (Vn/Vo) ' 7 ,V n <V 
\ MO , V n > Vq 

Here v n is the component of relative velocity along the line joining the grain centers, is the 
asymptotic coefficient at large velocities (typically 0.9 here), and vq = \/2ga, with g being the 
acceleration due to gravity and a the radius of the particle [17]. 

We examine three simulation geometries: "disks", a purely two-dimensional simulation 
with disks as used in [11]; "spheres in 2D" which are spherical particles in a three-dimensional 
(3D) simulation kept in a 2D plane by having reflecting front and back walls 2.002a apart 
and with reflecting left and right walls; and "spheres in 3D" which are spherical particles in 
a chute with periodic front and back walls 16 — 32a apart and with reflecting left and right 
walls. In all three types of simulations a sieve is located at the bottom as described above. 
For simulations it is most convenient to use units where g = 1, a = 1, and the grain mass 
M = 1. A typical simulation run is 10 4 — 10 5 time units, equivalent to 2 — 20 minutes of real 
time in a system made up of balls 3 mm in diameter. 

Steady-state density and velocity profiles for a system similar to that depicted in Fig.d a ) 
are shown in Fig.d D ) and (c). These profiles are averaged over the xz— cross-section of the 
32a x 32a chute. Three different regimes are evident. Grains initially placed at the top are 
little affected by collisions and accelerate uniformly with v y — —g, as shown in Fig. dd). 
This free-fall region ends abruptly just above the top of the dense column of grains seen in 
Fig. da). There is then a short fluid region extending about 25a between the free-fall and 
glassy region. The grains behave like a simple fluid in this region, with a parabolic profile of 
the v y velocity typical of Poiseuille flow (CPs in Fig. 2(b)). For nearly monodisperse systems, 
the density increases rapidly in this region until it reaches the random close-packed volume 
fraction of ~ 0.64. Except for a small region at the bottom, the remainder of the system has a 
density greater than 0.64 indicating partial crystallization. For systems of particles with 15% 
polydispersity, the density remains constant at 0.6 and is in a glass-like state. As a result we 
will now concentrate on 15% polydisperse systems as they are more uniform. 

Fig. db) shows the v y plug profile (B's) along with the corresponding shear stress profile 
in the glassy region of such a system. Clearly, the system supports a finite shear stress in the 
central region thus justifying calling this a glass. In steady state, the weight of the system, 
pg can be supported by either a pressure gradient d y P yy or by a gradient in the shear stress 
d x Pxy (pressures and stresses are measured as in [11]). A comparison of d y P yy and d x P xy (cf. 
Fig.sd e ) an d EJ shows that the shear stress, and hence the walls, carry most (~ 98%) of the 
weight. 

As jamming is approached, Donev et al. [4] showed that the pressure in a classical (con- 
servative) hard sphere glass approaches 

P = D{pk B T){l~p/ Pc )-\ (3) 

where D is the dimension, ks is Boltzmann's constant, and p c is the close-packed density. 
As can be seen in the inset in Fig.de), the diagonal components of the pressure tensor do 
approach this value as p approaches p c . However, there are small systematic deviations in the 
glassy phase, probably due to the dissipation and the fact that the velocities are not Maxwell 
distributed [18] as assumed in the derivation of Eq. 

It is reasonable to doubt the capability of a simple binary collision model to describe the 
dense glassy phase [19]. However, it is capable of reproducing a number of experimental results 
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Fig. 2 - Shear stress and velocity profiles for the 15% polydisperse 3D simulation. Data in the glassy 
region for h =50, 100, 150, 200 are plotted separately but fall on top of each other (B's). The D's in 
(b) show the velocity profile in the fluid region at h — 280 which gives a parabolic profile (solid line). 



related to this dense phase. For example, experiments have shown [9] that the fluctuating 
velocity i5v is related to the flow velocity v by a power-law relationship Sv oc v 2 / 3 . We 
observe the same relationship (Fig. |3J). Velocity fluctuations are not, however, isotropic. In 
a equilibrium fluid, fluctuations in the velocity are governed by equipartition and Sv 2 = 
(Vy) — (v y ) 2 is the same as Sv 2 and 5v 2 . This is very nearly the case in what we label the fluid 
region, but is not at all the case in the uncquilibrated free fall region or in the glassy region. 
Further measurements, such as that of the kurtosis of the velocity distribution, showed that 
the distribution of velocities is not Boltzmann-distributed. 

Considerable attention has also been focused on relating the force/impulse and collision 




Fig. 3 - Relationship between fluctuating and flow velocity in the glassy region. Data was averaged 
in directions normal to g for the 32x32 (*, Sv = f<k>| + Sv y + Sv^ 1 ^ 2 ) and 16x16 (A, Sv — Sv y ) 15% 
polydisperse systems. The fitted lines have slope of 2/3, in agreement with the experiments of [9]. v y 
is varied by varying the probability of reflection p from 0.01 to 0.995. 
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time distributions to jamming. A collision time is defined as the time between two consecutive 
collisions for a particle. In a simple fluid the collision times are distributed exponentially and 
the mean collision time (or "relaxation time") is proportional to the mean free path, the 
typical distance traveled between collisions. This is indeed what we observe in what we call 
the "fluid" region of the simulations, as shown in Fig.BJa). The distribution of collision times 
can be used to test the validity of many of the assumptions that go into standard kinetic 
theories used to describe granular gases and show us how they break down as we move into 
the dense glassy regions. As we will show, the distribution can also distinguish between a 
disordered glass and a crystalline packing of grains. 

The distributions of the collision times N(t) for binary collisions are shown in Fig. ^ For 
the glassy region of the 15% polydisperse systems (Fig.^Jb)), N(r) follows a power law, 

N(t) ~ t a (4) 

where a is 2.81 ±0.06 for distributions of grain-grain collisions in all the geometries, regardless 
of dimension. The power-law behavior of N(t) is independent of the specific choice made for 
the coefficient of restitution (as long as it is less than 1), and whether or not we use a constant 
or velocity-dependent coefficient of restitution as in Eq. We also examined the distribution 
of the distance particles travel between consecutive collisions. It has the same power-law 
behavior as the collision time distribution and we will therefore concentrate on the collision 
time distribution here. 

Experiments often refer to systems with less than 5% polydispersity as monodisperse. 
However, truly monodisperse systems eventually crystallize and give a power law with a = 4 
as shown in Fig. 0{c) (3D monodisperse systems did not crystallize over the entire system 
so the 3D data in Fig^Ic) is taken from a crystallized region). However, we find significant 
differences between even 1% polydisperse and truly monodisperse systems and systems with 
7.5% polydispersity yielded the same power laws as the 15% system. It appears that the 
polydispersity is a critical factor in these calculations mainly due to its affect on breaking up 
crystalization. 

The significance of the power-law distribution of collision times lies in the fact that it 
implies that collisions are not statistically independent events. If they were independent, it is 
straight-forward to show that the distribution of collision times would be exponential, or at 
least it would have an exponential tail [20] . This does not necessarily suggest a breakdown of 
the model based on binary collisions, suggested for other reasons in Ref. [19]. It does however 
suggest a breakdown in any kinetic theory arguments based on the Boltzmann equation which 
relies on the underlying assumption of molecular chaos and the statistical independence of 
collisions. Similarly, hydrodynamic descriptions which are based on a Chapman-Enskog ex- 
pansion of the Boltzmann distribution break down. This is due to the breakdown in the 
statistical independence of collisions, not because ternary collisions are required to describe 
the system. Thus, the collision time distribution can differentiate a granular fluid (exponen- 
tial), a granular crystal (power law a = 4), and a granular glass (power-law a = 2.8). 

Typically an exponential force distribution or impulse distribution is considered a signature 
of jamming. It has been suggested however [5-7] that an upturn at the lower end of the 
impulse distribution may provide more information. Impulse distributions from our spheres 
in 2D simulation, as shown in Fig. E{b), have a power law upturn at the low impulses and 
an exponential tail at the high impulses. The power law at the lower impulses appears to be 
more a signature of jamming as even the granular fluid phase exhibits an exponential tail. 

Our "spheres in 2D" simulation is a direct analog of the experiment described in [10], the 
main difference being that the experiment discharged the grains from a hopper at the bottom 
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Fig. 4 - (a) Unnormalized collision time distribution in the fluid region of a 15% polydisperse "3D 
spheres" simulation showing an exponential distribution, (b) Unnormalized collision time distribu- 
tions in the center of the glassy region for 15% polydisperse "disks" (dash-dotted line), "spheres in 
2D" (dotted line), and "3D spheres" (solid line) showing a power-law with exponent 2.85 as indicated 
by the thick gray line, (b) Similar to (a) but for monodisperse particles that have crystallized. The 
power-law in this case has exponent 4 as indicated by the thick gray line. Note, plots in (b) and (c) 
are log- log plots whereas (a) is a semi- log plot. All distributions are averaged over time in a region 
slightly larger than a grain size in the x and y directions and over the entire z direction. The precise 
choice of the location, other than it being clearly in the appropriate region, does not change the result. 

of the chute whereas we use a sieving process at the bottom. The experiment used a pressure 
transducer attached to the chute wall to detect grain-wall collisions, measure impulses, and 
times between collisions. However the experiment primarily sees an exponential distribution 
of impulses, although there is an upturn at small impulses for the slower flows. Further, for 
the distribution of grain- wall collision times they found a power-law with exponent a = 1.5. 
Based on the distribution of impulses, we explain the main discrepancy not as a surface versus 
bulk effect but as a result of experimental response time and sensitivity of the detector as 
follows. If we remove collisions with impulses in the power-law part of the distribution (i.e. 
the small impulse end that may be difficult to resolve experimentally) by putting in a cut off 




Fig. 5 - (a) Two-dimensional distribution of impulses and collision times (log scales) for a 15% 
"spheres in 2D" polydisperse simulation. The horizontal solid line in (a) indicates the cut-off in 
impulses, coinciding with the exponential tail of the ID impulse distribution shown in the semi-log plot 
(b), that results in the 1.5 collision time power-law observed in (c). Integrating the two-dimensional 
distribution in (a) over all collision times gives the ID histogram of impulses shown in (b) . Including 
only those impulses in the exponential tail of the impulse distribution (b), we obtain the 1.5 collision 
time power law as shown by the thick gray line in (c) . 
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at low impulses we also get the 1.5 collision time power law found in the experiment [10], as 
shown in Fig. [3] 

In conclusion, power laws with exponents independent of dimension were found for the 
distribution of collision times. Further, polydispersity and disorder are directly related to the 
power law exponent. The collision time distribution is an exponential for a granular fluid, 
a power-law with exponent a = 4 for a granular crystal, and a power-law with exponent 
a = 2.8 for a granular glass. It is interesting to note that the propagations of stress seems to 
show a similar ability to discern crystal and glass [3]. We successfully compared our collision 
time power law distribution with the experiment described in [10]. By subtracting out the 
collisions with very small impulses that would be difficult to measure experimentally, we found 
that the stronger collisions that make up the exponential tail of the impulse distribution 
have an associated collision time distribution with power law exponent 1.5. This highlights 
the importance of simulations that can give access to information not easily accessible in 
experiments. Measurement of the full collision-time distribution in experiments would require 
accessibility to the measurement of very low impulse collisions. In some systems this might 
be precluded by inelastic collapse and prolonged contacts. However, groups have observed 
peaks in the low end of the force distribution for static systems [5, 7] so it seems likely that 
a similar measurement of impulses in the analogous dynamic case could see the full spectrum 
of phenomenon we find in simulations. 
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